Measurement exploratory data analysis

DEVICES = ['icomox', 'icomox-fast', 'beaglebone', 'beaglebone-bvs', 'beaglebone-old', 'beaglebone-fan']

DEVICE = DEVICES[2]
T_WAVEFORM = 5
T_SEC = T_WAVEFORM
NFFT = 1024

from tabulate import tabulate
from IPython.display import Markdown, HTML
import matplotlib.pylab as plt
from typing import List, Tuple

import pandas as pd
import numpy as np
from scipy.signal import find_peaks
from scipy.signal import butter, lfilter

from tsfel.feature_extraction.features import fundamental_frequency

import os
import os.path
import sys
sys.path.append('../')
from vibrodiagnostics import mafaulda
from vibrodiagnostics import discovery
from vibrodiagnostics import datasets

pd.set_option("display.notebook_repr_html", False)

file_list = []

if DEVICE == 'icomox':
    Fs = 400
    DIR = '../../air_conditioning/shc3/'
    file_list = ['icomox_k5_low_power_accel.txt']
    file_list = [os.path.join(DIR, filename) for filename in file_list]
    DATASET = datasets.icomox_dataset(file_list)

elif DEVICE == 'icomox-fast':
    Fs = 50000
    DIR = '../../air_conditioning/shc3/'
    file_list = ['icomox_k5_accel.txt']
    file_list = [os.path.join(DIR, filename) for filename in file_list]
    DATASET = datasets.icomox_dataset(file_list)

elif DEVICE == 'beaglebone':
    Fs = 2500
    file_list = [
        'k3_1.tsv', 'k3_2.tsv', 'k3_3.tsv', 'k3_4.tsv',
        'k5_1.tsv', 'k5_2.tsv', 'k5_3.tsv', 'k5_4.tsv'
    ]
    DIR = '../../air_conditioning/shc3/'
    
    file_list = [os.path.join(DIR, filename) for filename in file_list]
    DATASET = datasets.beaglebone_dataset(file_list)

elif DEVICE == 'beaglebone-bvs':
    Fs = 2500
    file_list = [
        'bvs_1_hore.tsv', 'bvs_2_hore.tsv', 'bvs_3_motor.tsv', 'bvs_4_motor.tsv'
    ]
    DIR = '../../air_conditioning/bvs/'
    
    file_list = [os.path.join(DIR, filename) for filename in file_list]
    DATASET = datasets.beaglebone_dataset(file_list)

elif DEVICE == 'beaglebone-old':
    Fs = 1160     # 1160 Hz +/-150 Hz (1020 - 1340 Hz)
    file_list = [
        'excool_digitalis/high_pressure_pump_40percent__1.tsv',
        'excool_digitalis/high_pressure_pump_40percent__2.tsv',
        'excool_digitalis/high_pressure_pump_40percent_speed_up__3.tsv',
        'excool_digitalis/high_pressure_pump_80percent__4.tsv',
        'shc2/compressor_shc2_klima4_bad.tsv',
        'shc2/compressor_shc2_klima5_good.tsv',
        'vertiv_digitalis/compressor_top_1.tsv',
        'vertiv_digitalis/compressor_top_2.tsv',
        'vertiv_digitalis/compressor_top_3.tsv',
        'vertiv_digitalis/compressor_base.tsv',
        'vertiv_digitalis/compressor_side.tsv'
    ]
    DIR = '../../air_conditioning/'
    file_list = [os.path.join(DIR, filename) for filename in file_list]
    DATASET = beaglebone_dataset(file_list)

elif DEVICE == 'beaglebone-fan':
    Fs = 2500
    file_list = [
        'fan/speeds/1_still.tsv',
        'fan/speeds/2_still.tsv',
        'fan/speeds/3_still.tsv',
        'fan/speeds/1_up.tsv',
        'fan/speeds/2_up.tsv',
        'fan/speeds/3_up.tsv',
        'fan/speeds/1_down.tsv',
        'fan/speeds/2_down.tsv',
        'fan/speeds/3_down.tsv',
    ]
    DIR = '../../air_conditioning/'
    file_list = [os.path.join(DIR, filename) for filename in file_list]
    DATASET = beaglebone_dataset(file_list)

for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    ts.info()
    print()

k3_1.tsv

<class 'pandas.core.frame.DataFrame'>
Index: 51000 entries, 0.0 to 20.3996
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   x       51000 non-null  float64
 1   y       51000 non-null  float64
 2   z       51000 non-null  float64
dtypes: float64(3)
memory usage: 1.6 MB

k3_2.tsv

<class 'pandas.core.frame.DataFrame'>
Index: 51000 entries, 0.0 to 20.3996
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   x       51000 non-null  float64
 1   y       51000 non-null  float64
 2   z       51000 non-null  float64
dtypes: float64(3)
memory usage: 1.6 MB

k3_3.tsv

<class 'pandas.core.frame.DataFrame'>
Index: 51000 entries, 0.0 to 20.3996
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   x       51000 non-null  float64
 1   y       51000 non-null  float64
 2   z       51000 non-null  float64
dtypes: float64(3)
memory usage: 1.6 MB

k3_4.tsv

<class 'pandas.core.frame.DataFrame'>
Index: 51000 entries, 0.0 to 20.3996
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   x       51000 non-null  float64
 1   y       51000 non-null  float64
 2   z       51000 non-null  float64
dtypes: float64(3)
memory usage: 1.6 MB

k5_1.tsv

<class 'pandas.core.frame.DataFrame'>
Index: 51000 entries, 0.0 to 20.3996
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   x       51000 non-null  float64
 1   y       51000 non-null  float64
 2   z       51000 non-null  float64
dtypes: float64(3)
memory usage: 1.6 MB

k5_2.tsv

<class 'pandas.core.frame.DataFrame'>
Index: 51000 entries, 0.0 to 20.3996
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   x       51000 non-null  float64
 1   y       51000 non-null  float64
 2   z       51000 non-null  float64
dtypes: float64(3)
memory usage: 1.6 MB

k5_3.tsv

<class 'pandas.core.frame.DataFrame'>
Index: 51000 entries, 0.0 to 20.3996
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   x       51000 non-null  float64
 1   y       51000 non-null  float64
 2   z       51000 non-null  float64
dtypes: float64(3)
memory usage: 1.6 MB

k5_4.tsv

<class 'pandas.core.frame.DataFrame'>
Index: 51000 entries, 0.0 to 20.3996
Data columns (total 3 columns):
 #   Column  Non-Null Count  Dtype  
---  ------  --------------  -----  
 0   x       51000 non-null  float64
 1   y       51000 non-null  float64
 2   z       51000 non-null  float64
dtypes: float64(3)
memory usage: 1.6 MB

for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    display(tabulate(ts.describe(), headers='keys', tablefmt='html'))

k3_1.tsv

	x	y	z
count	51000	51000	51000
mean	2.30926e-15	1.14551e-15	3.23924e-16
std	1.37145	2.11749	0.512991
min	-2.46212	-5.23984	-1.22658
25%	-1.31252	-1.59941	-0.436227
50%	0.076596	-0.377948	-0.0530239
75%	1.39386	1.32252	0.47388
max	2.68717	4.69949	1.26424

k3_2.tsv

	x	y	z
count	51000	51000	51000
mean	-1.12739e-15	3.32283e-15	5.67347e-15
std	1.37462	1.90235	0.54732
min	-2.44795	-4.07667	-1.19419
25%	-1.32229	-1.4661	-0.475688
50%	0.0907731	-0.436238	0.00331546
75%	1.43198	1.07262	0.53022
max	2.48579	4.1143	1.36848

k3_3.tsv

	x	y	z
count	51000	51000	51000
mean	9.60069e-16	-1.59607e-15	3.9874e-16
std	1.40137	1.95904	0.55607
min	-2.62483	-4.38867	-1.25534
25%	-1.25967	-1.56255	-0.488938
50%	0.0336402	-0.269236	0.0140165
75%	1.4467	1.26358	0.540921
max	2.76396	4.16155	1.28338

k3_4.tsv

	x	y	z
count	51000	51000	51000
mean	-2.11937e-15	4.49606e-15	5.62652e-15
std	0.627253	0.8348	0.22793
min	-9.79898	-6.27935	-1.67302
25%	-0.0273018	-0.0283453	-0.0444018
50%	0.0205986	-0.0043951	-0.0204516
75%	0.0445488	0.0435053	0.0274488
max	7.75651	9.81518	1.72791

k5_1.tsv

	x	y	z
count	51000	51000	51000
mean	3.58434e-15	-2.53483e-15	6.39767e-16
std	1.49955	2.32511	0.231579
min	-3.22158	-5.17207	-0.665095
25%	-1.09002	-2.01065	-0.162141
50%	0.251196	-0.382032	-0.0184398
75%	1.2571	2.15669	0.149212
max	2.88572	4.95886	1.25092

k5_2.tsv

	x	y	z
count	51000	51000	51000
mean	-9.78738e-16	2.91908e-15	2.61535e-15
std	1.5101	2.32367	0.22692
min	-3.3225	-4.90644	-0.651701
25%	-1.09513	-2.05637	-0.148747
50%	0.270034	-0.451703	-0.0289962
75%	1.29989	2.13492	0.138655
max	2.97641	4.79339	1.19246

k5_3.tsv

	x	y	z
count	51000	51000	51000
mean	-5.92676e-15	-7.54011e-16	1.70753e-15
std	1.49944	2.33441	0.225151
min	-3.29453	-4.94197	-0.632079
25%	-1.09111	-2.044	-0.153075
50%	0.298001	-0.415387	-0.00937392
75%	1.25601	2.12333	0.134327
max	2.76487	4.90156	1.21209

k5_4.tsv

	x	y	z
count	51000	51000	51000
mean	8.0361e-15	8.27908e-15	4.04591e-16
std	1.40058	2.10618	0.258499
min	-4.15798	-5.68461	-1.33782
25%	-1.02051	-1.85257	-0.140307
50%	0.0572508	0.0155409	-0.0445065
75%	1.20686	1.85971	0.123145
max	5.87715	6.93715	3.76357

Time domain histogram

for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    axis = ts.columns
    ax = ts[axis].boxplot(figsize=(8, 3), grid=True)
    plt.show()

k3_1.tsv

k3_2.tsv

k3_3.tsv

k3_4.tsv

k5_1.tsv

k5_2.tsv

k5_3.tsv

k5_4.tsv

Time domain waveform

for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    axis = ts.columns
    
    ax = ts[axis].plot(figsize=(20, 8), grid=True, subplots=True)
    for i, axname in enumerate(axis):
        ax[i].set_xlabel('Time [s]')
        ax[i].set_ylabel(f'Amplitude ({axname}) [m/s^2]')
    plt.show()               # plt.savefig('waveform.png')

k3_1.tsv

k3_2.tsv

k3_3.tsv

k3_4.tsv

k5_1.tsv

k5_2.tsv

k5_3.tsv

k5_4.tsv

Time domain waveform zoom detail

for name, ts in DATASET:
    axis = ts.columns
    display(Markdown(f'**{name}**'))
    ax = (ts[axis].iloc[int(T_WAVEFORM*Fs):int(T_WAVEFORM*Fs)+Fs]
                  .plot(figsize=(20, 10), grid=True, subplots=True))
    
    for i, axname in enumerate(axis):
        ax[i].set_xlabel('Time [s]')
        ax[i].set_ylabel(f'Amplitude ({axname}) [m/s^2]')
        plt.show()      # plt.savefig('waveform_zoom.png')

k3_1.tsv

k3_2.tsv

k3_3.tsv

k3_4.tsv

k5_1.tsv

k5_2.tsv

k5_3.tsv

k5_4.tsv

def spectogram(x):
    fig, ax = plt.subplots(figsize=(15, 4))
    cmap = plt.get_cmap('inferno')
    pxx, freqs, t, im = plt.specgram(
        x, NFFT=NFFT, Fs=Fs,
        detrend='mean',
        mode='magnitude', scale='dB',
        cmap=cmap, vmin=-60
    )
    fig.colorbar(im, aspect=20, pad=0.04)
    ax.set_xlabel('Time [s]')
    ax.set_ylabel('Frequency [Hz]')
    mafaulda.resolution_calc(Fs, NFFT)
    return freqs, pxx


def window_idx(t):
    return (Fs * t) // NFFT + 1


def spectrum_slice(freqs, Pxx, t):
    fig, ax = plt.subplots(2, 1, figsize=(20, 8))
    n = window_idx(t)

    dB = 20 * np.log10(Pxx.T[n] / 0.000001)
    ax[0].plot(freqs, dB)      # 1 dB = 1 um/s^2
    ax[0].grid(True)
    ax[0].set_xlabel('Frequency [Hz]')
    ax[0].set_ylabel('Amplitude [dB]')
    
    ax[1].plot(freqs, Pxx.T[n])
    ax[1].grid(True)
    ax[1].set_xlabel('Frequency [Hz]')
    ax[1].set_ylabel('Amplitude [m/s^2]')
    return n


def get_max_frequency(freqs, Pxx, i):
    max_freq = freqs[np.argmax(Pxx.T[i])]
    return max_freq


def get_peaks(freqs, Pxx, i, top=5):
    amplitudes = Pxx.T[i]
    peaks, _ = find_peaks(amplitudes, distance=3)

    fundamental = get_max_frequency(freqs, Pxx, i)
    f_top = freqs[peaks[np.argsort(amplitudes[peaks])]][::-top]
    y_top = np.sort(amplitudes[peaks])[::-top]

    return pd.DataFrame({
        'f': f_top,
        'y': y_top,
        '1x': f_top / fundamental 
    })


def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
    b, a = butter(order, [lowcut, highcut], fs=fs, btype='band')
    y = lfilter(b, a, data)
    return y


def get_spectrograms(DATASET: List[pd.DataFrame], axis: str) -> list:
    spectrograms = []

    for name, ts in DATASET:
        base_freq = fundamental_frequency(ts[axis], Fs)
        display(Markdown(f'**{name}** *({axis.upper()} axis, Fundamental = {base_freq:.4f} Hz)*'))
        
        freqs, Pxx = spectogram(ts[axis])
        spectrograms.append((name, freqs, Pxx))
        plt.show()          # plt.savefig(f'x_axis_fft_{NFFT}.png')
    
    return spectrograms


def show_spectrogram_detail(spectrograms: list, axis: str, t: float):
    for name, freqs, Pxx in spectrograms:
        display(Markdown(f'**{name}** ({axis.upper()} axis @ {t}s)'))
        i_window = spectrum_slice(freqs, Pxx, t)
        plt.show()           #plt.savefig(f'x_axis_fft_{NFFT}_at_{T_SEC}s.png')


def show_mms_peaks(spectrograms: list, axis: str, t: float):
    for name, freqs, Pxx in spectrograms:
        display(Markdown(f'**{name}** ({axis.upper()} axis @ {t}s)'))
    
        i_window = window_idx(t)
        peaks = discovery.mms_peak_finder(Pxx.T[i_window])
        
        fig, ax = plt.subplots(1, 1, figsize=(15, 3))
        ax.grid(True)
        ax.plot(freqs, Pxx.T[i_window])
        ax.scatter(freqs[peaks], Pxx.T[i_window][peaks], marker='^', color='red')
        ax.set_xlabel('Frequency [Hz]')
        
        plt.show()


def show_harmonic_series(spectrograms: list, axis: str, t: float):
    # https://stackoverflow.com/questions/1982770/changing-the-color-of-an-axis
    for name, freqs, Pxx in spectrograms:
        display(Markdown(f'**{name}** ({axis.upper()} axis @ {t}s)'))
    
        i_window = window_idx(t)
        h_series = discovery.harmonic_series_detection(freqs, Pxx.T[i_window], Fs, NFFT)
    
        # Find best (sum of harmonics' amplitudes in the largest)
        max_harmonic_amp_idx = np.argmax([
            sum([h[1] for h in s]) / len(s)
            for s in h_series
        ])
        best_harmonic_series = pd.DataFrame(
            h_series[max_harmonic_amp_idx],
            columns=['Frequency [Hz]', 'Amplitude [m/s^2]']
        )
        best_harmonic_series.index += 1
        display(tabulate(best_harmonic_series, headers='keys', tablefmt='html'))
    
        # Plot found harmonic series
        fig, ax = plt.subplots(1, 8, figsize=(30, 4))
        for i in range(8):
            s = h_series[i+1]
            if i == max_harmonic_amp_idx:
                ax[i].xaxis.label.set_color('red')
    
            ax[i].plot(freqs, Pxx.T[i_window])
            ax[i].scatter([x[0] for x in s], [x[1] for x in s], marker='^', color='red')
            ax[i].set_xlabel('Frequency [Hz]')
    
        plt.show()

def show_spectra_largest_amplitudes(spectrograms: list, axis: str, t: float):
    for name, freqs, Pxx in spectrograms:
        display(Markdown(f'**{name}** ({axis.upper()} axis @ {t}s)'))

        i_window = window_idx(t)
        x_fundamental = get_max_frequency(freqs, Pxx, i_window)
        peaks = get_peaks(freqs, Pxx, i_window)
        
        display(Markdown(f'- *Fundamental frequency:* {x_fundamental} Hz'))
        display(tabulate(peaks.head(5), headers='keys', tablefmt='html'))

Spectrogram in X axis

x_spectra = get_spectrograms(DATASET, 'x')

k3_1.tsv (X axis, Fundamental = 47.7451 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

k3_2.tsv (X axis, Fundamental = 47.7451 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

k3_3.tsv (X axis, Fundamental = 47.7451 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

k3_4.tsv (X axis, Fundamental = 47.7941 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

k5_1.tsv (X axis, Fundamental = 48.1373 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

k5_2.tsv (X axis, Fundamental = 48.1373 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

k5_3.tsv (X axis, Fundamental = 48.1373 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

k5_4.tsv (X axis, Fundamental = 48.1373 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

Spectrogram detail in X axis

show_spectrogram_detail(x_spectra, 'x', T_SEC)

k3_1.tsv (X axis @ 5s)

k3_2.tsv (X axis @ 5s)

k3_3.tsv (X axis @ 5s)

k3_4.tsv (X axis @ 5s)

k5_1.tsv (X axis @ 5s)

k5_2.tsv (X axis @ 5s)

k5_3.tsv (X axis @ 5s)

k5_4.tsv (X axis @ 5s)

Peaks in frequency spectrum in X axis

• MMS peak finder algorithm

show_mms_peaks(x_spectra, 'x', T_SEC)

k3_1.tsv (X axis @ 5s)

k3_2.tsv (X axis @ 5s)

k3_3.tsv (X axis @ 5s)

k3_4.tsv (X axis @ 5s)

k5_1.tsv (X axis @ 5s)

k5_2.tsv (X axis @ 5s)

k5_3.tsv (X axis @ 5s)

k5_4.tsv (X axis @ 5s)

Harmonic series detection in X axis

# show_harmonic_series(x_spectra, 'x', T_SEC)

show_spectra_largest_amplitudes(x_spectra, 'x', T_SEC)

k3_1.tsv (X axis @ 5s)

• Fundamental frequency: 48.828125 Hz

	f	y	1x
0	48.8281	0.685933	1
1	285.645	0.180302	5.85
2	524.902	0.0620602	10.75
3	1052.25	0.0087056	21.55
4	717.773	0.00689525	14.7

k3_2.tsv (X axis @ 5s)

• Fundamental frequency: 48.828125 Hz

	f	y	1x
0	48.8281	0.703858	1
1	285.645	0.136529	5.85
2	524.902	0.0330748	10.75
3	954.59	0.00578816	19.55
4	14.6484	0.00447683	0.3

k3_3.tsv (X axis @ 5s)

• Fundamental frequency: 48.828125 Hz

	f	y	1x
0	48.8281	0.656088	1
1	285.645	0.149954	5.85
2	144.043	0.0424332	2.95
3	1164.55	0.0068831	23.85
4	639.648	0.00539692	13.1

k3_4.tsv (X axis @ 5s)

• Fundamental frequency: 4.8828125 Hz

	f	y	1x
0	4.88281	0.0768066	1
1	153.809	0.00184415	31.5
2	297.852	0.00152948	61
3	83.0078	0.00119425	17
4	1015.62	0.0010988	208

k5_1.tsv (X axis @ 5s)

• Fundamental frequency: 48.828125 Hz

	f	y	1x
0	48.8281	0.921647	1
1	288.086	0.102319	5.9
2	434.57	0.0270554	8.9
3	722.656	0.0109019	14.8
4	400.391	0.00575804	8.2

k5_2.tsv (X axis @ 5s)

• Fundamental frequency: 48.828125 Hz

	f	y	1x
0	48.8281	0.92099	1
1	288.086	0.10624	5.9
2	480.957	0.0224154	9.85
3	14.6484	0.00968361	0.3
4	627.441	0.00550303	12.85

k5_3.tsv (X axis @ 5s)

• Fundamental frequency: 48.828125 Hz

	f	y	1x
0	48.8281	0.936292	1
1	288.086	0.106108	5.9
2	241.699	0.0264785	4.95
3	295.41	0.0106427	6.05
4	1103.52	0.00580459	22.6

k5_4.tsv (X axis @ 5s)

• Fundamental frequency: 48.828125 Hz

	f	y	1x
0	48.8281	0.898246	1
1	288.086	0.0986729	5.9
2	480.957	0.027933	9.85
3	295.41	0.00986706	6.05
4	983.887	0.00442585	20.15

Spectrogram in Y axis

y_spectra = get_spectrograms(DATASET, 'y')

k3_1.tsv (Y axis, Fundamental = 47.7451 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

k3_2.tsv (Y axis, Fundamental = 47.7451 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

k3_3.tsv (Y axis, Fundamental = 47.7451 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

k3_4.tsv (Y axis, Fundamental = 47.3039 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

k5_1.tsv (Y axis, Fundamental = 48.1373 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

k5_2.tsv (Y axis, Fundamental = 48.1373 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

k5_3.tsv (Y axis, Fundamental = 48.1373 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

k5_4.tsv (Y axis, Fundamental = 48.1373 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

Spectrogram detail in Y axis

show_spectrogram_detail(y_spectra, 'y', T_SEC)

k3_1.tsv (Y axis @ 5s)

k3_2.tsv (Y axis @ 5s)

k3_3.tsv (Y axis @ 5s)

k3_4.tsv (Y axis @ 5s)

k5_1.tsv (Y axis @ 5s)

k5_2.tsv (Y axis @ 5s)

k5_3.tsv (Y axis @ 5s)

k5_4.tsv (Y axis @ 5s)

Peaks in frequency spectrum in Y axis

show_mms_peaks(y_spectra, 'y', T_SEC)

k3_1.tsv (Y axis @ 5s)

k3_2.tsv (Y axis @ 5s)

k3_3.tsv (Y axis @ 5s)

k3_4.tsv (Y axis @ 5s)

k5_1.tsv (Y axis @ 5s)

k5_2.tsv (Y axis @ 5s)

k5_3.tsv (Y axis @ 5s)

k5_4.tsv (Y axis @ 5s)

Harmonic series detection in Y axis

# show_harmonic_series(y_spectra, 'y', T_SEC)

show_spectra_largest_amplitudes(y_spectra, 'y', T_SEC)

k3_1.tsv (Y axis @ 5s)

• Fundamental frequency: 48.828125 Hz

	f	y	1x
0	48.8281	1.20267	1
1	478.516	0.198349	9.8
2	285.645	0.103564	5.85
3	573.73	0.0188074	11.75
4	622.559	0.0106902	12.75

k3_2.tsv (Y axis @ 5s)

• Fundamental frequency: 48.828125 Hz

	f	y	1x
0	48.8281	1.03286	1
1	144.043	0.138609	2.95
2	285.645	0.0364966	5.85
3	541.992	0.0092398	11.1
4	1240.23	0.00614032	25.4

k3_3.tsv (Y axis @ 5s)

• Fundamental frequency: 48.828125 Hz

	f	y	1x
0	48.8281	1.00388	1
1	478.516	0.12413	9.8
2	429.688	0.0633835	8.8
3	1164.55	0.0102918	23.85
4	1115.72	0.00770109	22.85

k3_4.tsv (Y axis @ 5s)

• Fundamental frequency: 9.765625 Hz

	f	y	1x
0	9.76562	0.0299383	1
1	21.9727	0.00412541	2.25
2	83.0078	0.00197565	8.5
3	244.141	0.00151702	25
4	227.051	0.0013613	23.25

k5_1.tsv (Y axis @ 5s)

• Fundamental frequency: 48.828125 Hz

	f	y	1x
0	48.8281	1.38965	1
1	385.742	0.147958	7.9
2	529.785	0.0615459	10.85
3	400.391	0.0299462	8.2
4	26.8555	0.0151213	0.55

k5_2.tsv (Y axis @ 5s)

• Fundamental frequency: 48.828125 Hz

	f	y	1x
0	48.8281	1.40826	1
1	434.57	0.179166	8.9
2	419.922	0.0564781	8.6
3	405.273	0.0231326	8.3
4	673.828	0.0110141	13.8

k5_3.tsv (Y axis @ 5s)

• Fundamental frequency: 48.828125 Hz

	f	y	1x
0	48.8281	1.42221	1
1	480.957	0.15385	9.85
2	578.613	0.0547047	11.85
3	402.832	0.0217529	8.25
4	1149.9	0.0119277	23.55

k5_4.tsv (Y axis @ 5s)

• Fundamental frequency: 48.828125 Hz

	f	y	1x
0	48.8281	1.3829	1
1	480.957	0.143644	9.85
2	578.613	0.0521242	11.85
3	983.887	0.0168376	20.15
4	14.6484	0.0122953	0.3

Spectrogram in Z axis

z_spectra = get_spectrograms(DATASET, 'z')

k3_1.tsv (Z axis, Fundamental = 47.7451 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

k3_2.tsv (Z axis, Fundamental = 47.7451 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

k3_3.tsv (Z axis, Fundamental = 47.7451 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

k3_4.tsv (Z axis, Fundamental = 47.7451 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

k5_1.tsv (Z axis, Fundamental = 0.0490 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

k5_2.tsv (Z axis, Fundamental = 48.1373 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

k5_3.tsv (Z axis, Fundamental = 48.1373 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

k5_4.tsv (Z axis, Fundamental = 48.1373 Hz)

Window size: 1024
Heinsenberg box
	Time step: 409.6 ms
	Frequency step: 2.44140625 Hz

Spectrogram detail in Z axis

show_spectrogram_detail(z_spectra, 'z', T_SEC)

k3_1.tsv (Z axis @ 5s)

k3_2.tsv (Z axis @ 5s)

k3_3.tsv (Z axis @ 5s)

k3_4.tsv (Z axis @ 5s)

k5_1.tsv (Z axis @ 5s)

k5_2.tsv (Z axis @ 5s)

k5_3.tsv (Z axis @ 5s)

k5_4.tsv (Z axis @ 5s)

Peaks in frequency spectrum in Z axis

show_mms_peaks(z_spectra, 'z', T_SEC)

k3_1.tsv (Z axis @ 5s)

k3_2.tsv (Z axis @ 5s)

k3_3.tsv (Z axis @ 5s)

k3_4.tsv (Z axis @ 5s)

k5_1.tsv (Z axis @ 5s)

k5_2.tsv (Z axis @ 5s)

k5_3.tsv (Z axis @ 5s)

k5_4.tsv (Z axis @ 5s)

Harmonic series detection in Z axis

# show_harmonic_series(z_spectra, 'z', T_SEC)

show_spectra_largest_amplitudes(z_spectra, 'z', T_SEC)

k3_1.tsv (Z axis @ 5s)

• Fundamental frequency: 48.828125 Hz

	f	y	1x
0	48.8281	0.280503	1
1	444.336	0.0312579	9.1
2	383.301	0.0217816	7.85
3	507.812	0.0143637	10.4
4	859.375	0.0133318	17.6

k3_2.tsv (Z axis @ 5s)

• Fundamental frequency: 48.828125 Hz

	f	y	1x
0	48.8281	0.32505	1
1	429.688	0.0253058	8.8
2	144.043	0.0191175	2.95
3	954.59	0.0163078	19.55
4	1049.8	0.0117508	21.5

k3_3.tsv (Z axis @ 5s)

• Fundamental frequency: 48.828125 Hz

	f	y	1x
0	48.8281	0.323556	1
1	334.473	0.0279629	6.85
2	668.945	0.0219133	13.7
3	456.543	0.0152852	9.35
4	922.852	0.010976	18.9

k3_4.tsv (Z axis @ 5s)

• Fundamental frequency: 4.8828125 Hz

	f	y	1x
0	4.88281	0.0167918	1
1	61.0352	0.00451767	12.5
2	90.332	0.00312322	18.5
3	324.707	0.00191779	66.5
4	424.805	0.00173085	87

k5_1.tsv (Z axis @ 5s)

• Fundamental frequency: 48.828125 Hz

	f	y	1x
0	48.8281	0.0713598	1
1	434.57	0.0264955	8.9
2	192.871	0.0225645	3.95
3	34.1797	0.0192975	0.7
4	1010.74	0.0169037	20.7

k5_2.tsv (Z axis @ 5s)

• Fundamental frequency: 48.828125 Hz

	f	y	1x
0	48.8281	0.0682579	1
1	192.871	0.0260875	3.95
2	336.914	0.0217514	6.9
3	144.043	0.0200225	2.95
4	31.7383	0.0125301	0.65

k5_3.tsv (Z axis @ 5s)

• Fundamental frequency: 48.828125 Hz

	f	y	1x
0	48.8281	0.0773615	1
1	144.043	0.0257784	2.95
2	7.32422	0.0197851	0.15
3	288.086	0.0169507	5.9
4	422.363	0.0140544	8.65

k5_4.tsv (Z axis @ 5s)

• Fundamental frequency: 48.828125 Hz

	f	y	1x
0	48.8281	0.0907712	1
1	144.043	0.0249169	2.95
2	625	0.0199538	12.8
3	964.355	0.0172125	19.75
4	14.6484	0.0139624	0.3

Histogram

axis = ['x', 'y', 'z']
for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    ts[axis].hist(figsize=(10, 5), grid=True, bins=50)
    plt.show()

k3_1.tsv

k3_2.tsv

k3_3.tsv

k3_4.tsv

k5_1.tsv

k5_2.tsv

k5_3.tsv

k5_4.tsv

axis = ['x', 'y', 'z']
for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    ts[axis].boxplot(figsize=(10, 5))
    plt.show()

k3_1.tsv

k3_2.tsv

k3_3.tsv

k3_4.tsv

k5_1.tsv

k5_2.tsv

k5_3.tsv

k5_4.tsv

Orbitals of all cross sections

for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    fig, ax = plt.subplots(1, 3, figsize=(20, 4))

    for i, col in enumerate([('x', 'y'), ('x', 'z'), ('y', 'z')]):
        ax[i].scatter(ts[col[0]], ts[col[1]], s=1)
        ax[i].grid(True)
        ax[i].set_xlabel(col[0].upper())
        ax[i].set_ylabel(col[1].upper())
        ax[i].grid(True)
    plt.show()       # plt.savefig('orbitals.png')

k3_1.tsv

k3_2.tsv

k3_3.tsv

k3_4.tsv

k5_1.tsv

k5_2.tsv

k5_3.tsv

k5_4.tsv

Orbitals of 1x harmonic frequency

x_spectra_by_name = {spec[0]: spec for spec in x_spectra}
y_spectra_by_name = {spec[0]: spec for spec in y_spectra}
z_spectra_by_name = {spec[0]: spec for spec in z_spectra}
t = 5
space = 5

for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    fig, ax = plt.subplots(1, 3, figsize=(20, 4))

    name, freqs, Pxx = x_spectra_by_name[name]
    x_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))
    name, freqs, Pxx = y_spectra_by_name[name]
    y_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))

    name, freqs, Pxx = z_spectra_by_name[name]
    z_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))

    try:
        ts['x_1x'] = butter_bandpass_filter(ts['x'], x_fundamental - space, x_fundamental + space, Fs)
        ts['y_1x'] = butter_bandpass_filter(ts['y'], y_fundamental - space, y_fundamental + space, Fs)
        ts['z_1x'] = butter_bandpass_filter(ts['z'], z_fundamental - space, z_fundamental + space, Fs)
    except ValueError:
        continue
    
    for i, col in enumerate([('x_1x', 'y_1x'), ('x_1x', 'z_1x'), ('y_1x', 'z_1x')]):
        ax[i].scatter(ts[col[0]], ts[col[1]], s=1)
        ax[i].grid(True)
        ax[i].set_xlabel(col[0].upper())
        ax[i].set_ylabel(col[1].upper())
        ax[i].grid(True)
    
    plt.show()       # plt.savefig('orbitals_1x.png')

k3_1.tsv

k3_2.tsv

k3_3.tsv

k3_4.tsv

k5_1.tsv

k5_2.tsv

k5_3.tsv

k5_4.tsv

t = 5
space = 8

for name, ts in DATASET:
    display(Markdown(f'**{name}**'))

    name, freqs, Pxx = x_spectra_by_name[name]
    x_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))
    name, freqs, Pxx = y_spectra_by_name[name]
    y_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))

    name, freqs, Pxx = z_spectra_by_name[name]
    z_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))

    try:
        x = butter_bandpass_filter(ts['x'], x_fundamental - space, x_fundamental + space, Fs)
        y = butter_bandpass_filter(ts['y'], y_fundamental - space, y_fundamental + space, Fs)
        z = butter_bandpass_filter(ts['z'], z_fundamental - space, z_fundamental + space, Fs)
    except ValueError:
        continue

    ax = plt.figure().add_subplot(projection='3d')
    ax.scatter(x, y, z, zdir='z', s=1, color='navy')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_zlabel('z')
    ax.set_xlim(-2, 2)
    ax.set_ylim(-2, 2)
    ax.set_zlim(-2, 2)
    ax.zaxis.labelpad = -0.7
    plt.show()

k3_1.tsv

k3_2.tsv

k3_3.tsv

k3_4.tsv

k5_1.tsv

k5_2.tsv

k5_3.tsv

k5_4.tsv